/**
 * mLibMat33.js
 * Created: Oct 2013
 * Source link: http://code.google.com/p/code-depository/source/browse/
 *
 * By Suresh Nalluri,
 * License: The MIT License (MIT)
 *          See LICENSE details
 *
 *
 *	Todo:
 *
 *
 * Requires:
 *	mLibFactory.js
 *	mLibIncludes.js
 *	mLibVec3.js
 *	mLibQuaternion.js,
 *
 *
 * 3x3 Matrix Class
 * A matrix object is composed of a one dimensional array of 9 elements.  The matrix is sorted in
 * column major order, which imples that element [i,j] is the (i+3*j)th entry in the array.
 *
 *
 * note this constructor pattern allows us to create objects in 2 ways:
 * 1. var o = mLib.Mat33( 0 or more params );      //OR
 * 2. var o = new mLib.Mat33( 0 or more params );
 *
 * NOTE: var o = [new] mLIB.Mat33(<no params>); // will return default identity matrix *
 *
 * added benefit of this pattern is that we get a copy constructor by default!!
 * eg. var newMtx = mLIB.Mat33(existingMtx);
 * so no need to have an explicit copy method.
 *
 * Usage:
 *		var l = [new] mLIB.Mat33( mat | quat |  <noargs> ]);
 * 
 *-------------------------------------------------------
 **/
if (typeof mLIB == 'undefined') function mLIB() {};

( function(module) {

    var eqZero = module.Includes.isZero,
        sqrt = module.Includes.sqrt,
        sin = module.Includes.sin,
        cos = module.Includes.cos,
        acos = module.Includes.acos,
        atan2 = module.Includes.atan2,
        PI = module.Includes.PI,
        ZERO = module.Includes.ZERO,
        axisIdx = ['x', 'y', 'z'];
    
	module.Mat33 = Factory.obj( {
        init: function() {
            var A = this.mat = new Array(9);
            
            if (arguments[0].isClass(module.Mat33) ) {
                // param is a matrix - copy constructor
                var src = arguments[0].mat;
                
                A[0] = src[0]; A[3] = src[3]; A[6] = src[6];
                A[1] = src[1]; A[4] = src[4]; A[7] = src[7];
                A[2] = src[2]; A[5] = src[5]; A[8] = src[8];
                
            } else if (arguments[0].isClass(module.Quat) {
                // param is a quaternion - axis-angle construction
                this.rotationByQuaternion(arguments[0]);
                
            } else if (arguments[0] == 0) {
                // return the zero matrix
                A[0] = A[3] = A[6] = 0;
                A[1] = A[4] = A[7] = 0;
                A[2] = A[5] = A[8] = 0;                

            } else {
                // identity matrix
                A[0] = 1; A[3] = 0; A[6] = 0;
                A[1] = 0; A[4] = 1; A[7] = 0;
                A[2] = 0; A[5] = 0; A[8] = 1;
                  
            }
        }, //init()
        
        // synonym for the '=' operand, eg. mtx.equals(other) - copies other into mtx.
        equals: function(m) {
            if ( this === m ) return;
            
            var A = this.mat,
                B = M.mat;
            
            A[0] = B[0]; A[3] = B[3]; A[6] = B[6];
            A[1] = B[1]; A[4] = B[4]; A[7] = B[7];
            A[2] = B[2]; A[5] = B[5]; A[8] = B[8];
            return this;
        },

        CopyOf: function(m) {
            // ie. clone
            return module.Mat33( m );
        },
        
        // clean self.. set near-zero elements to zero.
        clean: function() {
            var A = this.mat;
            
            for (var i = 0; i < 9; i++) {
                eqZero(A[i]) && (A[i] = 0.0);
            }
        },
        
        Clean: function() {
            return module.Mat33(this).zero();
        },
        
        // set self to identity matrix
        identity: function() {
            var A = this.mat;
            A[0] = 1; A[3] = 0; A[6] = 0;
            A[1] = 0; A[4] = 1; A[7] = 0;
            A[2] = 0; A[5] = 0; A[8] = 1;
            return this;
        },
       
        // comparison operations
        isEqual: function(m) {
            if ( this === m )return true;  // since it's the same matrix object!!
            
            var A = this.mat,
                B = m.mat;
                
            for (var i = 0; i < 9; i++) {
                if ( !eqZero(A[i] - B[i]) return false;
            }
            return true;
        },
        
        isNotEqual: function(m) {
            var A = this.mat,
                B = m.mat;
                
            for (var i = 0; i < 9; i++) {
                if ( !eqZero(A[i] - B[i]) return true;
            }
            return false;
       },
       
        isZero: function() {
            var A = this.mat;                

            for (var i = 0; i < 9; i++) {
                if ( !eqZero(A[i]) return false;
            }
            return true;
        },
        
        isIdentity: function() {
            var A = this.mat;
            
            return (
                eqZero(1.0-A[0]) && eqZero(A[3]) && eqZero(A[6]) &&
                eqZero(A[1]) && eqZero(1.0-A[4]) && eqZero(A[7]) &&
                eqZero(A[2]) && eqZero(A[5]) && eqZero(1.0-A[8])
            )
        
        },
        
        // transpose self
        transpose: function() {
            var tmp,
                A = this.mat;
            
            tmp = A[3];
            A[3] = A[1];
            A[1] = tmp;
            
            tmp = A[6];
            A[6] = A[2];
            A[2] = tmp;
            
            tmp = A[7];
            A[7] = A[5];
            A[5] = tmp;
            
            return this;
        },

        // return the transpose of self in a new Mat33
        Transpose: function() {
            return new module.Mat33(this).transpose();
        }
        
        // set self to inverse
        inverse: function() {
            var A = this.mat,
                coFactor0 = A[4]*A[8] - A[5]*A[7],
                coFactor3 = A[2]*A[7] - A[1]*A[8],
                coFactor6 = A[1]*A[5] - A[2]*A[4],
                det = A[0]*coFactor0 + A[3]*coFactor3 + A[6]*coFactor6;  // the determinant
            
            if (eqZero(det)) return false;
            
            // inverse = 1/det * [Adjoint] - Cramer's method
            //
            var recip = 1 / det,
                M = this.toArray(); // make a copy
            
            A[0] = recip*coFactor0;
            A[1] = recip*coFactor3;
            A[2] = recip*coFactor6;
            
            A[3] = recip*(M[5]*M[6] - M[3]*M[8]);
            A[4] = recip*(M[0]*M[8] - M[2]*M[6]);
            A[5] = recip*(M[2]*M[3] - M[0]*M[5]);
            
            A[6] = recip*(M[3]*M[7] - M[4]*M[6]);
            A[7] = recip*(M[1]*M[6] - M[0]*M[7]);
            A[8] = recip*(M[0]*M[4] - M[1]*A[3]);
        },
        
        Inverse: function() {
            return module.Mat33(this).inverse();
        },
        
        determinant: function() {
            var A = this.mat;
            
            return A[0]*(A[4]*A[8] - A[5]*A[7]) +
                   A[3]*(A[2]*A[7] - A[1]*A[8]) +
                   A[6]*(A[1]*A[5] - A[2]*A[4]);
        },
        
        Adjoint: function() {
            var M = this.mat,
                adjoint = module.Mat33(),
                B = out.mat;

            //get the cofators and transpose them
            B[0] = M[4]*M[8] - M[5]*M[7]
            B[1] = M[2]*M[7] - M[1]*M[8]
            B[2] = M[1]*M[5] - M[2]*M[4]
            
            B[3] = M[5]*M[6] - M[3]*M[8]
            B[4] = M[0]*M[8] - M[2]*M[6]
            B[5] = M[2]*M[3] - M[0]*M[5]
            
            B[6] = M[3]*M[7] - M[4]*M[6]
            B[7] = M[1]*M[6] - M[0]*M[7]
            B[8] = M[0]*M[4] - M[1]*M[3]
            
            return adjoint;
        },
        
        trace: function() {
            return (this.mat[0] + this.mat[4] + this.mat[8]); // sum of the diagonal elements
        },
        
        multiplyByScalar: function(f) {
            var A = this.mat;
            A[0] *= f; A[3] *= f; A[6] *= f;
            A[1] *= f; A[4] *= f; A[7] *= f; 
            A[2] *= f; A[5] *= f; A[8] *= f;
            return this;
        },
        
        MultiplyByScalar: function(f) {
            return module.Mat33(this).multiplyByScalar(f);
        },
        
        plus: function(m) {
            var A = this.mat,
                B = m.mat;
                
            for (var i =0; i< 9; i++) {
                A[i] += B[i];
            }
            return this;
        },
        Plus: function(m) {
            var out = module.Mat33(this), // make a copy
                A = m.mat,
                B = out.mat;
            
            for (var i =0; i< 9; i++) {
                B[i] += A[i];
            }
            return out;

        },
        
        minus: function(m) {
            var A = this.mat,
                B = m.mat;
                
            for (var i =0; i< 9; i++) {
                A[i] -= B[i];
            }
            return this;
        },

        Minus: function(m) {
            return module.Mat33(this).plus(m);
        },

        negation: function() {
            var A = this.mat;
            
            for (var i =0; i< 9; i++) {
                A[i] = -A[i];
            }
            return this;
        
        },
        
        Negation: function() {
            return module.Mat33(this).negation();
        },
        
        multiply: function(m) {
            var A = this.toArray(),  // cheap than creating a copy of this!!
                B = m.mat,
                C = this.mat;
        
            C[0] = A[0]*B[0] + A[3]*B[1] + A[6]*B[2];
            C[3] = A[0]*B[3] + A[3]*B[4] + A[6]*B[5];
            C[6] = A[0]*B[6] + A[3]*B[7] + A[6]*B[8];
            
            C[1] = A[1]*B[0] + A[4]*B[1] + A[7]*B[2];
            C[4] = A[1]*B[3] + A[4]*B[4] + A[7]*B[5];
            C[7] = A[1]*B[6] + A[4]*B[7] + A[7]*B[8];
            
            C[2] = A[2]*B[0] + A[5]*B[1] + A[8]*B[2];
            C[5] = A[2]*B[3] + A[5]*B[4] + A[8]*B[5];
            C[8] = A[2]*B[6] + A[5]*B[7] + A[8]*B[8];
                    
            return this;
        },
        
        Multiply: function(m) {
            return module.Mat33(this).multiply(m);
        },
        
        // Matrix - Column vector multiply ie c = [A]b
        multiplyColumn: function(colVec) {
            var A = this.mat,
                b = colVec.toArray();       // this may not be cheap !!
            
            return module.Vec3( 
                A[0]*b[0] + A[3]*b[1] + A[6]*b[2],
                A[1]*b[0] + A[4]*b[1] + A[7]*b[2],
                A[2]*b[0] + A[5]*b[1] + A[8]*b[2]
            );
        },    
        
        // Row vector - Matrix multiply  ie. c = b[A]
        multiplyRow: function(rowVec) {
            var A = this.mat,
                b = rowVec.toArray();
            
            return module.Vec3(
                b[0]*A[0] + b[1]*A[1] + b[2]*A[2],
                b[0]*A[3] + b[1]*A[4] + b[2]*A[5],
                b[0]*A[6] + b[1]*A[7] + b[2]*A[8]
            );
        },
        
        // set the matrix based on the given quaternion
        rotationByQuaternion: function(Q /*quaternion*/) {
            // assumtion: Q is a unit quaternion
            var xs = Q.x+Q.x,
                ys = Q.y+Q.y,
                zs = Q.z+Q.z,
                wx = Q.w*xs,
                wy = Q.w*ys,
                wz = Q.w*zs,
                xx = Q.x*xs,
                xy = Q.x*ys,
                xz = Q.x*zs,
                yy = Q.y*ys,
                yz = Q.y*zs,
                zz = Q.z*zs,
                A = this.mat;

            A[0] = 1.0 - (yy + zz);
            A[3] = xy - wz;
            A[6] = xz + wy;
            
            A[1] = xy + wz;
            A[4] = 1.0 - (xx + zz);
            A[7] = yz - wx;
            
            A[2] = xz - wy;
            A[5] = yz + wx;
            A[8] = 1.0 - (xx + yy);

            return this;
        
        },
        
        // set the matrix to a rotation matrix - given the rotation around each axis individually
        // ie rotation based on euler angles
        //
        rotationByEuler: function(xAngle, yAngle, zAngle) {
            // if the rotations by axis are represented by the matrices: Rx, Ry, Rz, then
            // the rotation = RxRyRz
            var cx = cos(xAngle), sx = sin(xAngle),
                cy = cos(yAngle), sy = sin(yAngle),
                cz = cos(zAngle), sz = sin(zAngle),
                A = this.Mat;
            
            A[0] =  cy*cz;            A[3] = -cy*sz;            A[6] = sy;
            A[1] =  sx*sy*cz + cx*sz; A[4] = -sx*sy*sz + cx*cz; A[7] = -sx*cy;
            A[2] = -cx*sy*cz + sx*sz; A[5] =  cx*sy*sz + sx*cz; A[8] = cx*cy;
            
            return this;
        },
        
        // set rotation around an arbitrary axis.. by axis we mean vector
        // around which we rotate
        rotationByAxisAngle: function(axis, angle) {
            var c = cos(angle), s = sin(angle), p = 1.0 - c,
                n = axis.Normalise(),       // normalised vector in the axis direction
                px = p * n.x, py = p * n.y, pz = p * n.z,
                sx = s * n.x, sy = s * n.y, sz = s * n.z,
                pxy = px * n.y, pyz = py * n.z, pxz = px * n.z,
                A = this.mat;
                
            A[0] = px*n.x + c; A[3] = pxy - sz;   A[6] = pxz + sy;
            A[1] = pxy + sz;   A[4] = py*n.y + c; A[7] = pyz - sx;
            A[2] = pxz - sy;   A[5] = pyz + sx;   A[8] = pz*n.z + c;

            return this;
        },
        
        // set self as rotation matrix, based on angle in radians around the x-axis
        rotationX: function(theta) {
            var sinTheta = sin(theta),
                A = this.mat;
            
            // +-            -+
            // | 1.0 0.0 0.0  |
            // | 0.0 cos -sin |
            // | 0.0 sin cos  |
            // +-            -+
            
            A[0] = 1.0;
            A[1] = A[2] = A[3] = A[6] =0.0;
            A[4] = A[8] = cos(theta); A[5]= sinTheta;
            A[7] = -sinTheta;
            
            return this;
        },
        
         // set self as rotation matrix, based on angle in radians around the y-axis
        rotationY: function(theta) {
            var sinTheta = sin(theta),
                A = this.mat;

            // +-            -+
            // | cos  0.0 sin |
            // | 0.0  1.0 0.0 |
            // | -sin 0.0 cos |
            // +-            -+

            A[0] = A[8] = cos(theta);
            A[2] = -sinTheta; A[6] = sinTheta;
            A[1] = A[3] = A[5] = A[7] 0.0;
            A[4] = 1.0;
            
            return this;
        },
         // set self as rotation matrix, based on angle in rdians around the z-axis
        rotationZ: function(radians) {
            var sinTheta = sin(theta),
                A = this.mat;

            // +-            -+
            // | cos -sin 0.0 |
            // | sin  cos 0.0 |
            // | 0.0  0.0 1.0 |
            // +-            -+
            
            A[0] = A[4] = cos(theta);
            A[1] = sinTheta;
            A[3] = -sinTheta;
            A[2] = A[5] = A[6] = A[7] = 0.0;
            A[8] = 1.0;
            
            return this;
        },
        
        // get one set of possible z-y-x fixed angles that will generate this matrix
        // Assumes that this is the fixed angles (z-x-y) rotation matrix:
        //  +-                                         -+
        //  |  cy*cz             -cy*sz           sy    |   where cy = Cos(yTheta)
        //  |  sx*sy*cz+cx*sz    -sx*sy*sz+cz*cz -sx*cy |
        //  | -cx*sy*cz+sx*sz     cx*sy*sz+sx*cz  cx*cy |
        //  +-                                         -+
        //
        // This is also the Euler angles (roll-pitch-yaw or x-y-z) matrix.
        //
        // note that what we return is one possible solution out of many due to
        // using inverse trig functions, which can result in multiple results
        //
        getFixedAngles: function(zyx) {
            var sx, sy, sz,
                cx, cy, cz,
                A = this.mat;
            
            sy = A[6];
            cy = sqrt(1.0 - sy*sy);
            if ( eqZero(cy) ) {
                // we have gimbal lock problem: x and z axes are aligned, which implies
                // we can't tell which is rotation around x and which is rotation around
                // y.  So let's assume that rotation around z is zero
                sz = 0.0;
                cz = 1.0;
                sx = A[5];
                cx = A[4];
                
            } else {
                // no gimbal lock
                var f = 1.0 / cy;
                sx = -A[7]*f;
                cx = A[8]*f;
                sz = -A[3]*f;
                cz = A[0]*f;
                
            }
            
            // return the angles in the given array param.
            // note use of atan2 as it returns angles within the range: [-PI, PI]
            zyx[0] = atan2(sz, cz); // the zTheta
            zyx[1] = atan2(sy, cy); // the yTheta
            zyx[2] = atan2(sz, cz); // the xTheta
        
            return this;
        },
        
        // Gets one possible axis-angle pair that will generate this matrix
        // Assumes that this is the rotation matrix around an arbitrary axis r:
        //  +-                                     -+
        //  | p*x^2+c       p*x*y-s*z   p*x*z+s*y   |   p = 1 - cos(theta)
        //  | p*x*y+s*z     p*y^2+c     p*y*z-s*x   |   c = cos(theta), s = sin(theta)
        //  | p*x*z-s*y     o*y*z+s*x   t*z^2+c     |   rAxis = (x, y, z)
        //  +-                                     -+
        //  based on http://www.geometrictools.com/Documentation/RotationIssues.pdf
        //
        getAxisAngle: function(rAxis) {
            // we return the axis in the given Vec3 type parameter rAxis
            // and the angle theta in the return statement.
        
            var A = this.mat,
                trace = A[0] + A[4] + A[8],     // = 2Cos(theta) + 1
                theta = acos( (trace - 1.0) * 0.5 );
            
            if (eqZero(theta)) {
                // theta is zero, so we can return any old axis!
                // in our case we return the x Axis
                rAxis.x = module.Vec3.xAxis.x;
                rAxis.y = module.Vec3.xAxis.y;
                rAxis.z = module.Vec3.xAxis.z;
                return 0.0;
            }

            if (theta < PI-ZERO) {
                // our standard case, theta lies in range (0, PI)
                rAxis.x = A[5] - A[7];
                rAxis.y = A[6] - A[2];
                rAxis.z = A[1] - A[3];
                return theta;
            }
            
            // theta == PI
            var i = 0;
            
            if ( A[4] > A[0] )
                i = 1;
            if ( A[8] > A[i + 3*i] )
                i = 2;
            
            var j = (i+1)%3,
                k = (j+1)%3,
                s = sqrt( A[i + 3*i] - A[j + 3*j] - A[k + 3*k] + 1.0 ),
                f = 1.0 / s;
            
            rAxis[axisIdx[i]] = 0.5 * s;
            rAxis[axisIdx[j]] = (A[i + 3*j])*f;
            rAxis[axisIdx[k]] = (A[k + 3*i])*f;

            return theta;
        },
        
        // set self as scaling matrix based on a vector of scaling factors
        // ie. set the diagnonal elements of the matrix to the vector and the rest as 0's 
        scaling: function(sf) {
            var A = this.mat;
            
            A[0] = sf.x;
            A[4] = sf.y;
            A[8] = sf.z;
            A[1] = A[2] = A[3] = A[5] = A[6] = A[7] = 0.0;
            return this;
        },
    
        /**
         * Accessor methods.. setters and getters
         **/

        // set matrix row by row
        setRows: function( inRow1, inRow2, inRow3 ) {
            // another form of constructor!  the row's are of type Vec3
            // note we overwrite the existing mtx
            var A = this.mat;                

            A[0] = inRow1.x; A[3] = inRow1.y; A[6] = inRow1.z;
            A[1] = inRow2.x; A[4] = inRow2.y; A[7] = inRow2.z;
            A[2] = inRow3.x; A[5] = inRow3.y; A[8] = inRow3.z;
            return this;
        },
        
        getRows: function(out) {
            // out is an object, which gets its keys set
            // better to get an object parameter since, the calling app can reuse
            // rather than we create each time.. expensive!!
            var A = this.mat;                

            out.row1 = module.Vec3( A[0], A[3], A[6] );
            out.row2 = module.Vec3( A[1], A[4], A[7] );
            out.row3 = module.Vec3( A[2], A[5], A[8] );
            return this;  // for good measure
        }
        
        // return the jRow of the matrix as a vector. jRow is 0, 1 or 2
        getRow: function(jRow) {
            if (jRow > 2) return null;
            
            var A = this.mat;                
            return module.Vec3(A[jRow], A[3+jRow], A[6+jRow]);
        },
        
        //set matrix column by column
        setColumns: function( inCol1, inCol2, inCol3 ) {
            // another form of constructor!  the col's are of type Vec3
            // note we overwrite the existing mtx
            var A = this.mat;                

            A[0] = inCol1.x; A[3] = inCol2.x; A[6] = inCol3.x;
            A[1] = inCol1.y; A[4] = inCol2.y; A[7] = inCol3.y;
            A[2] = inCol1.z; A[5] = inCol2.z; A[8] = inCol3.z;
            return this;
        },
        
        getColumns: function(out) {
            var A = this.mat;                
            out.col1 = module.Vec3( A[0], A[1], A[2] );
            out.col2 = module.Vec3( A[3], A[4], A[5] );
            out.col2 = module.Vec3( A[6], A[7], A[8] );
            return this;
        },
        
        getColumn: function(iCol) {
            if (iCol) > 2) return null;
            
            var A = this.mat,
                i = 3*iCol;
                
            return module.Vec3( A[i], A[i+1], A[i+3] );
        },
        
        getElement: function(iRow, jCol) {
            return this.mat[iRow + 3*jCol];
        },


        $reader: function() {
            return ( function(iRow, jCol) {
                return this.mat[iRow + 3*jCol];
            }).bind(this);
        },

        // return the raw matrix array for fast access (direct) in other functions
        $rawArray: function() {
            return this.mat;
        },
 
        toArray: function() {
            var A = this.mat,
                arr = new Array(9);
            
            for (var i = 0; i< 9; i++) {
                arr[i] = A[i];
            }
            return arr;
        },
        
        log: function() {
            var A = this.mat;

            return [
                '[ ',
                [A[0], ', ', A[3], ', ', A[6], ', '].join(''),
                [A[1], ', ', A[4], ', ', A[7], ', '].join(''),
                [A[2], ', ', A[5], ', ', A[8]].join(''),
                ' ]'
            ].join('\n');
        }

    } ); // module.Mat33 = Factory.obj()

    // friend methods ---------------------------------------------
    //

    var friends = {
        // multiply m by scalar and return m
        multipyByScalar: function(m, f) {
            var B = m.mat;
                
            B[0] *= f; B[3] *= f; B[6] *= f;
            B[1] *= f; B[4] *= f; B[7] *= f; 
            B[2] *= f; B[5] *= f; B[8] *= f;
            
            return m;
        },
        
        // transpose m and return it
        transpose: function(m) {
            var A = m.mat,
                tmp

            tmp = A[3];
            A[3] = A[1];
            A[1] = tmp;
            
            tmp = A[6];
            A[6] = A[2];
            A[2] = tmp;
            
            tmp = A[7];
            A[7] = A[5];
            A[5] = tmp;
            
            return m;
        },
        
        // negate m and return it
        negation: function(m) {
            var A = m.mat;
            
            for (var i =0; i< 9; i++) {
                A[i] = -A[i];
            }
            return m;
        },

        // return a skew symetric matrix based on a 3d vector. 
        // we'll probably never use this since calculating cross products is done more efficiently in the 
        // vector library. we add it here for completion to show that it could be replaced using a skew
        // symmetric matrix instead
        //
        getSkewSymmetricMat: function(v) {
            // this matrix is useful for calculating the cross product vxw
            // ie. vxw = Sw, where S is the matrix
            //          0 Vy
            //      
            var S = Mat(0),     // get the zero matrix
                m = S.mat;
                
            m[1] = v.z; 
            m[2] = -v.y;
            m[3] = -v.z;
            m[5] = v.x;
            m[6] = v.y;
            m[7] = -v.x;
            
            return S;
        },

        toArray: function(M, arr) {
            var A = M.mat;
            
            for (var i = 0; i < 9; i++) {
                arr[i] = A[i];
            }
            return arr;
        },

        $reader: function(m) {
            return function(iRow, jCol) {
                return m.mat[iRow + 3*jCol];
            }
        },

        // return the raw matrix array for fast access (direct) in other functions
        $rawArray: function(m) {
            return m.mat;
        },
         
    }; // friends{}

    for (var f in friends) {
        friends.hasOwnProperty(f) && (module.Mat33[f] = friends[f]);
    };

} )(mLIB);